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The vapor-liquid critical behavior of intrinsically asymmetric fluids is studied in finite systems of 
linear dimensions, L, focusing on periodic boundary conditions, as appropriate for simulations. The 
recently propounded "complete" thermodynamic [L — > oo) scaling theory incorporating pressure 
mixing in the scaling fields as well as corrections to scaling [arXiv:cond-mat/0212145] , is extended 
to finite L, initially in a grand canonical representation. The theory allows for a Yang- Yang anomaly 
in which, when L — > oo, the second temperature derivative, (d 2 u a /dT 2 ), of the chemical potential 
along the phase boundary, p a (T), diverges when T — > T c — . The finite-size behavior of various special 
critical loci in the temperature-density or (T, p) plane, in particular, the fc-inflection susceptibility 
loci and the Q-maximal loci — derived from Ql{T, (p)l) = {m 2 )\ / '(m 4 ) l where m = p — (p)l - 
is carefully elucidated and shown to be of value in estimating T c and p c . Concrete illustrations are 
presented for the hard-core square-well fluid and for the restricted primitive model electrolyte in- 
cluding an estimate of the correlation exponent v that confirms Ising-type character. The treatment 
is extended to the canonical representation where further complications appear. 



I. INTRODUCTION AND OVERVIEW 

True phase transitions arise in statistical mechanics 
only in the thermodynamic limit in which the volume of 
a system, V = L d (in d dimensions), and the number of 
particles in the system, N, go to infinity, while the den- 
sity p = N/V remains finite. In this limit, to be denoted 
for brevity simply by L — > oo, the free energy and other 
quantities may exhibit singularities at a phase boundary 
or critical point as functions of the temperature or other 
thermodynamic fields. However, for finite systems as, 
in particular, realized in computer simulations, the free 
energy becomes analytic everywhere in the temperature 
and in other fields such as the chemical potential, p., and 
the pressure, p. Thus thermodynamic quantities that 
vary discontinuously or diverge in the thermodynamic 
limit become rounded when L is finite. 

Computer simulations have been useful in quantify- 
ing and gaining insights into phase transitions in vari- 
ous systems. Nevertheless, to obtain precise, sharp re- 
sults from simulations — inevitably performed on finite- 
systems — one must perform appropriate extrapolations 
on the size, L, of the simulation "box." Crucial ques- 
tions then arise: How should one best estimate critical 
points from the finite-size data? And, especially: How 
can one reliably ascertain the critical universality class 
of particular model systems? 

To study the statistical mechanics of finite systems, 
one must at the start address two basic issues, namely, 
the overall geometry of the system and the specific nature 
of the boundary conditions. Here we will have in mind 
general c£-dimensional systems with periodic boundary 
conditions imposed on "rectangular" boxes of dimensions 
L\ x L2 x • ■ ■ x Ld = V = L d in which the ratios L^jL 
remains fixed (typically at 1) when L 00. Of course, 
this geometry combined with periodic boundary condi- 
tions has been used extensively in computer simulations 



for studies of the bulk properties of fluids. 

In the case of critical phenomena in systems with a 
well defined axis of symmetry in some thermodynamic 
plane, notably model magnetic materials and analogous 
lattice gases [1], in which the critical density is trivially 
known and the variation with (T — T c ) is of primary in- 
terest, the long-established theory of finite-size scaling 
[2,3] and its subsequent developments [4-6], has provided 
effective answers to many questions of how to extrapo- 
late data for finite systems. However, two new issues 
that demand further consideration have recently come to 
the fore. These are, first, the desire to obtain precise, 
unbiased answers for the universal critical behavior of 
"complex" and, especially, asymmetric fluid systems — 
in which, in particular, both the critical temperature T c , 
and the critical density p c , must be accurately estimated 
[7] — and, second, the realization that the existence of a 
so-called Yang- Yang anomaly [8,9] — in which the chem- 
ical potential p. a {T) on the vapor- liquid phase boundary 

exhibits a divergent curvature when T — > T c requires 

a significant elaboration [8,10,11] of earlier formulations 
of bulk, thermodynamic scaling for fluids [12,13]. 

The appropriately extended, "complete" scaling for- 
mulation for bulk properties (i.e., in the thermodynamic 
limit) that is needed to encompass a Yang- Yang anomaly 
[8] has recently been carefully expounded and investi- 
gated in some detail: first in Part I of this article [10], to 
be denoted here as I, and, more fully, in the thesis [11] 
of the first author, which will be referred to here as K. It 
proves necessary to "mix" the pressure, p, into the linear 
(and nonlinear) scaling fields [8]. To be explicit, let us, 
following I, introduce the dimensionless deviations from 
the (bulk) critical point (p c ,T c ,/i c ) via 
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component fluid must, in general, take the forms 

p = p - k t - ZoA) (I- 2 ) 
t = t-hp-jip, (1.3) 
h = p - kit - j 2 p, (1.4) 

in which the quadratic and higher order terms have been 
dropped (see I). The crucial new feature (going beyond 
the previously accepted analyses: see [12,13]) is the pres- 
ence of the, in general, nonzero dimensionless pressure- 
mixing coefficients ji and j 2 ■ when these vanish the ear- 
lier formulations are satisfactory. 

In terms of the (nonlinear) scaling fields the general 
scaling hypothesis of I asserts that the thermodynamics 
near criticality can be described, at least asymptotically, 

by 

y(\ 2 - a p, At, X A h; \- 9i u 4 , \- 9s u 5 , ■ • •) = 0, (1.5) 

where A is a free, positive scaling parameter. The expo- 
nents a (for the specific heat) and A (for the ordering 
field h), are related to the other standard critical expo- 
nents via 

A = 2-a-/? = /? + 7 = /3<5, (1.6) 

while 84 = 9 and #5 are the positive leading even and 
odd correction-to-scaling exponents for the correspond- 
ing irrelevant scaling fields, m{p,T,p) and u 5 (p,T,p). 
One then discovers [8,1, K] that the scaling form (1.5) 
implies: (a) the existence of a Yang- Yang anomaly in 
which (d 2 p a I dT 2 ) diverges as ~ J 2 / 1 * | " when t — > and 
(b) a leading singular term varying as ~ j 2 \t\ 2lS in the co- 
existence curve diameter, that dominates the previously 
known term ~ (Zi + .;/i)|t| 1 ~ Q since, e.g., [3 = 0.326 and 
a = O.IO9 for d = 3 Ising-type criticality. Further new, 
singular terms of similar character appear in other ther- 
modynamic properties: see I. 

The task addressed here, in Sec. II, is to systematically 
extend the general formulation for bulk scaling, as cm- 
bodied in (1.1)-(1.5), to finite systems characterized by 
a (single) finite length scale, L. According to the gen- 
eral principles of finite-size scaling, by which all lengths 
should, in the critical region, be scaled by the correlation 
length, £(T) ~ we may anticipate that, in effect, 

the scaling parameter A in (1.5) may, in a grand canon- 
ical setting, be replaced by L x ' v . Let us also note that 
when, as for real fluids, hyperscaling is valid (see Sec. 
II. A) we have 

dv = 2-a. (1.7) 

It has, however, been pointed out [14] that the scaling 
fields, p, t, h, • ■ ■, themselves may, in a finite system, 
gain an explicit dependence on the size L. Thus finite- 
size effects in a system confined by hard walls might well 
be dominated by 1/L contributions [5,15]. This issue, 
which is by no means definitively settled in general, is 



considered briefly in Sec. II with the conclusion that for 
the case of periodic boundary conditions, which is our 
main concern here, one should anticipate additive terms 
in (1.2)-(1.4); specifically, then, we will [setting l n = 1; 
see 1(3.22)] adopt the scaling field 

pip, T,p;L)=p-k t-p- s /L d + • • • , (1.8) 

and likewise, with new coefficients si and s 2l for 
t(p,T,p;L) and h(p,T, p; L), with d > 2. (Note that 
the coefficients s , Si, and s 2 carry dimensions of L d .) 
Fortunately, it then transpires that these L-dependcnt 
contributions do not enter the leading behavior of the 
quantities of principal interest, such as the k- and Q-loci 
in the (T, p) plane: see below. 

Specific predictions for the finite-size variation of basic 
densities and susceptibilities are presented in Sec. II. C. 
The variation with L of the chemical potential, p, at 
the bulk critical temperature and density is examined in 
Sec. II. D: the answer provides a route to uncovering the 
presence of L-dependent terms in the scaling fields as in 
(1.8). 

Now, as mentioned, an important application of finite- 
size scaling theory is to analyze numerical data obtained 
from simulations on finite systems, and, thereby, to gain 
knowledge of the critical properties of the bulk system. 
Major efforts have been devoted to estimating critical pa- 
rameters such as (T c ,p c ), and to confirming universality 
classes. As regards the estimation of T c and p c , most 
studies have focused on calculating the coexistence curve 
in the (p, T) plane and then fitting the data with some 
suitably chosen formula in which T c and p c appear. 

However, simulations of a system in its two-phase re- 
gion may require prohibitively long times or special, more 
elaborate computational techniques to equilibrate the 
two coexisting phases owing to the free-energy barrier 
that grows rapidly as T decreases and L increases. More- 
over, since the correlation length, £(T,p), becomes large 
and eventually diverges when the critical region is ap- 
proached, finite-size effects smear out the vapor and liq- 
uid states near T c and blur their distinction thereby seri- 
ously hampering the reliable determination of the coex- 
istence curve. Finally, field mixing (even in the absence 
of pressure mixing) distorts the shape of the diameter, 
etc. Consequently, naively fitting coexistence curve data 
may yield quite poor values for T c and p c . 

To meet these latter challenges, Bruce and Wilding 
some time ago [16,17] proposed a rather convenient and 
effective finite-size scaling method for estimating T c and 
pc, that, in particular, incorporates fi and t mixing into 
the scaling fields t and h (although pressure mixing is 
not included). Their method, which has proved quite 
popular, is based on the hypothesis that fluid critical- 
ity belongs to the Ising universality class (or, more gen- 
erally, to some well studied universality class for which 
certain detailed critical properties are well established 
numerically). On that basis their method matches dis- 
tribution functions of density and energy fluctuations ob- 
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served in simulations to the (presumed available) limiting 
fixed point distributions as obtained a priori from simu- 
lations of simpler models (known to be of Ising or other 
character). In this way, following extrapolation on L, 
they estimate critical parameters. However, significant 
questions remain: What should be done when a priori 
knowledge of the (suspected or, possibly, quite new) crit- 
ical behavior of the system of interest is not available? 
And; How should one proceed if the effects of pressure 
mixing may not be negligible? [18] 

In light of these serious issues, an important aim of 
our studies has been to develop unbiased finite-size scal- 
ing methods for estimating T c and p c without the need 
for such strong assumptions and extensive a priori knowl- 
edge. For this purpose, as previously reported [7,19,10], 
various special loci have been introduced that, in the 
thermodynamic limit, spring from the critical point in 
the density-temperature or other thermodynamic plane. 
The bulk scaling behavior of these critical loci was de- 
rived within the complete, scaling theory in I (and also 
studied there within classical mean- field theory) . Among 
these loci, the fc-loci — defined via the points of isother- 
mal maxima of x^ = XI ' P k m the (p, T) plane, where 
X = p 2 k^TKT is the isothermal susceptibility — have 
already been used in simulations to estimate the criti- 
cal points of the hard-core square- well (HCSW) fluid [7], 
and of the restricted primitive model (RPM) electrolyte 
[19]. It is a goal of the present article to analyze the 
behavior of these fc-loci in systems of finite size: explicit 
expressions for p( k \T;L), the fc-loci, in the (p,T) plane 
are obtained in Sec. III. Not surprisingly, one finds that 
the density p = p^(T c ;L) evaluated on a fc-locus at T c 
(where we suppose that T c has been estimated reliably in 
some other way) approaches the critical density p c when 
L — > oo: But in what manner? 

We show in Sec. III. A that there is a leading devi- 
ation of magnitude L~ 2/3 / 1/ followed by a term of or- 
der L~^~°^l v : however, the amplitude of the leading 
contribution vanishes when k takes an "optimal" value 
fcopt = 372.^. In this result IZ^ is the (dimensionless) 
strength of the Yang- Yang anomaly as defined in Ref. [8] 
and in I. Sec. III. E. Extrapolating data for the densities 
p( k \T c ; L) to the thermodynamic limit can thus provide 
unbiased (bulk) estimates of the critical density. In Sec. 
III.B we re-apply this approach to the HCSW fluid us- 
ing what we believe is an improved estimate for T c : see 
below. Our new estimate for p c agrees well, within the 
uncertainties, with the previous result [7]. As indicated, 
this method for estimating p c has also been successfully 
applied to the RPM electrolyte [19]. 

Evidently, however, in locating p c by this route, one 
first needs a good estimate of T c . For fluids with rel- 
atively weak asymmetry like the hard-core square-well 
model, it was found [7] that the extrema in density of 
the fc-loci themselves provide fairly good estimators for 
T c that may be extrapolated in L. However, the whole 
critical region of the RPM is extremely asymmetric, in 
part, so it seems, because of the remarkably low value, 



p* = p c a? ~ 0.08 [19], of the reduced critical density 
(where a is the hard-core diameter). As a result, esti- 
mators for T c based on the available fc-loci prove rather 
misleading: indeed, the fc-loci for "near-optimal" values 
of k are observed to vary nonmonotonically in p — prob- 
ably as a result of competition between the two leading 
contributions, Ap ~ (k - k opt )/L 2fi / v and l/L^-^/", 
mentioned above. To overcome this serious obstacle to 
progress, Luijten, Fisher and Panagiotopoulos [19] intro- 
duced the Q-loci which they defined by points of isother- 
mal maxima in the (p, T) plane of the inverse Binder 
parameter [20] 

Ql{T;(p)l)= K t ^T with m = p-(p) L , (1.9) 

where (-)l denotes a grand-canonical ensemble average 
in the finite system. 

Now when L — > oo anywhere in the one-phase region 
one has Ql{T; p) — > \ [20], where for brevity, we have re- 
placed the argument (p) l in Ql by p. On the other hand, 
at criticality, Ql{T c ; p c ) approaches a universal value Q c 
that is close to 0.6326 for (d = 3)-dimensional Ising sys- 
tems in a cubic box with periodic boundary conditions 
[21-23]. For finite systems at fixed T near criticality, 
however, one finds that Ql exhibits rounded maxima 
that serve to provide well-defined loci, Pq(T; L) [19]. The 
behavior of these Q-loci for large L is derived explicitly 
within the full finite-size scaling theory in Sec. IV. A. One 
might note that determining the Q-loci involves calcu- 
lation of the fourth density moment and of its density 
derivative (i.e., the fifth moment) so that the analysis re- 
quires some care! By the same token, in order to obtain 
the Q-loci reliably via simulations, data of high quality 
are needed. As for the fc-loci, one may define Q^-loci by 
points of isothermal maxima in the (p, T) plane of a mod- 
ified Q parameter, namely, = Q L /p k . The behavior 
of these loci is presented in Sec. IV. B: we find that the 
density, Pq\t c ;L), evaluated at T c on these loci varies 

in leading order as L~ 2/3 / u with, as in the fc-loci, a sub- 
sequent i^f 1- ")/^ term. However, the amplitude of the 
leading contribution now vanishes when k = —9TZ^ 7 in 
contrast to fc Q pt = 37£ M for the fc-loci; thus the "optimal" 
value of k for the Q( fe )-loci has the opposite sign. 

Following Binder's original approach for symmetric 
systems [20], Luijten et al. [19] examined plots of 

Ql(T) = Q L (T;p Q (T;L)), (1.10) 

i.e., Ql evaluated on the Q-loci pq(T,L). For the RPM 
they observed that the successive self-intersections as 
L increased, say T®(L), converged rather rapidly to a 
precisely defined value, Too — which thus served as a 
good estimate for T c . At the same time, they surpris- 
ingly found that the values of Q® at the intersection 
points approached a limit which could be identified as a 
(surprisingly precise) estimate of the universal value Q c . 
Thereby they established convincingly that the RPM (at 
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least within the £ = 5 level of discretization they studied 
[19]) belongs to the short-range Ising universality class 
- despite the long-range Coulomb interactions in the 
model. We show here that the approach of the estima- 
tors, T®(L), derived from the Ql(T) plots to T c obeys 
a l/L( 1+0 y» law, while the difference Q%(T C (L)) - Q c 
varies as L~ 6 I V followed by a j|L _2 ^/ I/ term (see Sec. 
V.B). Note that these results are independent of asym- 
metry or pressure mixing (in leading order). 

In Sec. V.A we develop the theory for this approach 
and apply it to re-estimate T c for the HCSW model [7]. 
The new estimate is about 0.06% higher than the ear- 
lier value [7]; but that leads to no significant changes in 
the main conclusions reached previously: in particular, 
as noted above, the previous estimate for p c remains un- 
changed (within the uncertainties). 

On the other hand, in Sec. IV. C we consider the behav- 
ior of Ql{T; (p)l) for large L in the two-phase region be- 
neath T c . (See also Rovere, Heermann and Binder [24].) 
We exhibit plots for the HCSW fluid and RPM that illus- 
trate some striking features (and we correct a misleading 
expression given in [7] for the behavior of Ql(T; p^) with 
Pl = {p)l when L — > oo below T c ). In Sec. IV. D we go 
on to discuss the explicit scaling form for two minima of 
Ql{T; pl) that, when T < T c , approach the two sides of 
the coexistence curve rather rapidly as L — > oo: see Figs. 
8 and 9, below. It turns out that these considerations 
lead to a novel and apparently very effective and system- 
atic method of estimating the limiting coexistence curve 
width and diameter, namely, 

A Poo (T) = p + (T) - p_(T), (1.11) 
p d (T) = ±[p + (T) + p_(T)], (1.12) 

where p+{T) = pu q (T) and p~(T) = p va , P (T) denote the 
true, bulk liquid and vapor densities, respectively. This 
method, which yields precise results surprisingly close to 
T c , has been applied to the HCSW and RPM models; 
however, the details, which entail using the simulation 
data to generate a scaling function for the minima as 
T — > T c — , will be expounded elsewhere [25]. 

The universality class of a particular system can be 
identified or checked and confirmed by determining crit- 
ical exponents, a, (3, etc. In Sec. V.C we analyze further 
a method for estimating the correlation- length exponent, 
v [7]. This method has been applied to the HCSW fluid 
[7] and, more recently, reported for the RPM electrolyte 
[19]. A thermodynamic quantity for a finite system, say 
Pl(T), evaluated on some suitable locus, say p = p c , may 
exhibit a maximum at T = T C P (L) which can be regarded 
as an effective finite-size critical temperature. According 
to finite-size scaling one expects T C P (L) to approach the 
true critical temperature T c asymptotically as L~ x l v . We 
confirm that this conclusion survives pressure-mixing (for 
suitable loci) and, by way of an application, show that 
by examining a rather wide range of properties Pl(T) 
for the RPM one can identify those for which the desired 
maxima approach T c from above. This is important in 



practice because simulations above criticality are signif- 
icantly less hampered by problems of full equilibration 
than those at or below T c where two distinct putative 
phases coexist, and "alternate" in the simulation box. 
Consequently, sufficiently precise calculations of T^{L) 
are relatively easy which, in turn, provides a suitable ba- 
sis for robust extrapolation. In this way, we show that 
one can estimate the exponent v fairly accurately. For 
the RPM electrolyte (at the £ = 5 level of discretization) 
we find v = 0.63 ± 0.03 [19] which supports the conclu- 
sion that the model belongs to the (d — 3)-dimensional 
Ising universality class [19]. 

Both for gaining insight into experiments, in which 
the density p is most often a controlled variable, and, 
likewise, for simulations in which the particle number, 
N, is fixed, it is valuable to study the finite-size scaling 
behavior of near-critical fluids in a canonical or (p, T) 
representation. The bulk canonical free energy density 
f(p,T) = limL^ F N (V,T)/V, where F N (V,T) is the 
Helmholtz free energy, has a leading asymptotic scaling 
behavior near criticality of the form 

f(p, T) « f Q (p, T) + A\t\-V-^X ± (m/\tf), (1.13) 

in which /o (p, T) is a smooth (generally analytical) back- 
ground part of the free energy while m = (p — p c )/p c . 
However, this simple scaling form does not incorporate 
any mixing in the scaling fields. We may anticipate 
that upon incorporating the mixing of the scaling fields, 
the leading scaling behavior remains unchanged but with 
some modifications of the scaling variables m and t. But 
what should be expected precisely! That may well af- 
fect the behavior of the corrections on various loci [26]. 
And what scaling form should one obtain if, in particu- 
lar, pressure-mixing is introduced? In Sec. VI we derive 
explicit canonical scaling forms from the complete scaling 
formulation in the grand canonical representation. This 
is carried out first for the thermodynamic limit: then 
our finite-size results are applied to obtain corresponding 
canonical expressions. In Sec. VLB we discuss the defini- 
tion of finite-size canonical critical points and elucidate 
their behavior as illustrated by results for the HCSW 
fluid and the RPM electrolyte [7,19]. 

Finally, Sec. VII summarizes the article briefly. 

II. FULL FINITE-SIZE SCALING 
FORMULATION 

Here we extend to finite systems near bulk critical 
points the complete scaling theory that incorporates pres- 
sure mixing [10]. 

A. Scaling functions and hyperuniversality 

To extend the bulk scaling ansatz (1.5) to a finite 
V=L d system we first replace p, h, t; 114, U5, • • • by corre- 
sponding finite-size nonlinear scaling fields p(p,T, p; L), 
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Uj{p,T, /z; L), ■■■ of the form (1.8), etc, and 
choose an arbitrary fixed reference length, say /*. Set- 
ting A = (L/U) 1 /" in (1.5) then leads to the general 
hypothesis 

■--)=»• 

(2.1) 

which we expect to be at least asymptotically valid for 
L/l* — > oo as p, i and h — > 0. 

Let us now restrict attention to dimensionalities d less 
than the upper critical dimensionality o!> (= 4 for normal 
fluid criticality) . Then the hyperuniversality exponent re- 
lation, supported by renormalization group (RG) theory 
(for a fixed point without dangerous irrelevant variables 
[27]) dictates (2 — a)/v = d [see (1.7)] and we may solve 
(2.1) for p to obtain 

p c p(p, T, fi; L) = L- d Y(x L ,y L ; y L4 , • • •), (2.2) 

where we have introduced the dimensionless scaled vari- 
ables 

x L =D L tL 1 l v , y L = U L hL A /», y Lk = U Lk L- e «/» 

(fc = 4,5,---). (2.3) 

Here Dl, Ul, and Ulu oc u k are nonuniversal metrical 
factors, of dimensions I* ', Z„ A / u ^ l^l" ^ respectively, 
which depend on the system under study. 

By construction (note the factor p c > 0) the scaling 
function Y(x, y; y&, ■ ■ •) is dimensionless [28]. However, 
the hyperuniversality scaling hypothesis [29] (supported 
by various exact calculations [29-31], simulations [32] and 
RG theory [21]) tells us that Y(x,y;y4, • • •) is a univer- 
sal function of its (appropriately normalized) arguments. 
Note, however, that Y(x, y; 2/4, • • •) must depend on the 
geometry of the finite system and on the boundary condi- 
tions imposed; but it will not depend on any microscopic 
details beyond those that determine the bulk universality 
class of the relevant critical point. Furthermore, Y must 
be even under change of sign of the odd scaling variables, 

V -y, 2/5 -2/5, • • •• 

The bulk limit may be obtained formally by setting 
L = \j\DjJ\ v and letting L — > 00 (when it drops out 
of the nonlinear scaling fields p, t, • • •). This yields the 
scaling form 1(2.3), namely, 

p = Q\t\ 2 - a W ± (y; yi ,y 5 ,---), (2.4) 

with the identification Q = Z?i| 2-Q /p c , which is, thus, 
a dimensionless nonuniversal amplitude, while 

W±(y;y 4 ,---) = Y(±l,y;y 4 ,---), (2.5) 

is universal with the amplitudes in 1(2.1) and 1(2.2) re- 
lated by U= U L /\D L \ A , U k =U Lk \D L \ e - (fc = 4, 5, • • •)• 



In contrast to the bulk scaling function, the finite-size 
function, Y{x, y; 2/4, • • •) must be analytic in the vicinity 
of the origin since all critical singularities will be rounded 
in a finite system. Following I we may thus expand for 
large L in powers of the irrelevant variables to obtain 

Y(x L ,y L ; ■■■) = Y°(x L ,y L ) + £ Y"(x L , y L )yW , (2.6) 

where, as in I, the multi-index, k, is defined by k = (4), 
(5), • • •, (4,4), (4,5), • • •, (4,4,4), • • •, while y[^">™] means 
VhiVLj • • • Vhn- The underlying symmetry of the scaling 
function, Y(x L ,y L ; ■ ■ •), that is evidenced by exact results 
and RG theory, then requires 

Y K (x L ,-y L ) = ±Y K (x L ,y L ), (2.7) 

for k even or odd in the sense of 1(2.7). Thence we have 
the expansions 

Y"(x Ll y L ) = Y£ + Yf x L + Y^x\ + Y&y\ + • • • , 

= y L (Y£ + Y£x L + Y«x\ + Y^y\ + •••), 

(2.8) 

for At, even and odd, respectively, where the expansion 
coefficients Yu are universal numbers. 

For our present purposes the leading approximation 

Y « Y°(x L ,y L ) + y c L4 Y^ (x L ,y L ) + y c L5 Y^(x L , y L ), 

(2.9) 

in which Ula and Ul5 in the definitions of yLA and 3/1,5 
have been replaced by their critical-point values, will am- 
ply suffice. 



B. Finite-size corrections to the scaling fields 

In this section we discuss in a little more detail the 
question of finite-size corrections to the scaling fields that 
was touched on in the Introduction. This issue seems to 
have been first raised in Ref. [14] but to have escaped 
much more extensive or systematic discussion. Here we 
consider only a o!-dimensional hypercube with periodic 
boundary conditions. 

A field-theoretic RG approach to finite-size scaling was 
initiated by Brezin [33]. Later, with Zinn-Justin [21] sys- 
tematic calculations of the scaling functions were pre- 
sented using both d = 4 — e and d = 2 + e expansions. 
In particular, the shift of T c that enters the scaling vari- 
able, t, of the universal scaling functions was computed: 
see Ref. [21] Eqns. (3.20) and (3.32). Indeed, t as cal- 
culated in Eq. (3.21) of Ref. [21] contains finite-size cor- 
rections that, in leading order, vary as L~ 2 . A similar 
form for t was obtained by Korutcheva and Tonchev [34] 
for a finite system with long-range interactions decay- 
ing as i / / r d + 2 - 2o ' j (j _ > o+. Recently, Chen and Dohm 
[35] calculated the finite-size free-energy density of an 
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O(n) (p 4 field theory confined in a hypercube with peri- 
odic boundary conditions: they used a sharp cutoff in k 
space and obtained a nonuniversal L~ 2 contribution that 
dominated a universal scaling part that varied as L~ d . 

On the other hand, Jasnow and coworkers [31,36] con- 
cluded via RG theory that the system size L does not 
enter in the formation of the scaling fields: see, especially 
Ref. [31] Sec. III. Likewise Zinn-Justin [6, page 778] ar- 
gues that: "The crucial observation which explains finite- 
size scaling is that the renormalization theory which leads 
to RG equations is completely insensitive to finite size 
effects since renormalizations are entirely due to short 
distance singularities. As a consequence RG equations 
are not modified. ■ •" . Nevertheless, in our assessment 
it remains uncertain whether or not, even in the simplest 
case of periodic boundary conditions, the system size af- 
fects the scaling fields. While further careful analyses 
may settle the issue convincingly, we feel justified in al- 
lowing for an L" d leading contribution in all the scaling 
fields — as embodied in (1.8); however, it seems safe to 
assume that d > 2. As mentioned in the Introduction, we 
then find in most cases that these corrections are less im- 
portant, when L becomes large, than those arising from 
field mixing and the leading irrelevant variables. 

C. Some basic thermodynamic properties 

The generalized number and entropy "scaling" densi- 
ties, p and s, introduced in I play a significant role also 
in analyzing finite systems: they are defined by 

p = (dp/dh) h S = (dp/dt)- h . (2.10) 

From (2.2) and (2.3), we obtain [28] 

Pc p = U L L-f""(dyY), Pc S = D L L-^/-(d x Y), 

(2.11) 

where, here and below, we adopt the notations (d x Y) = 
(dY/dxL)y L , etc. 

Now recall the definitions 1(2.14) of the "true" reduced 
number and entropy densities, namely, 

p=£-=(®P) x _ S _ ( d i>\ (212) 
~ Pc \dpj t 7 p c k B 

Following I(2.16)-(2.19) these may be expressed in terms 
of the generalized, scaling densities. Thus we find 

p = l + (2q + l n )p + (n + 2l m )p + (v + l n 3 )t 
+ (1 - hh)p - {h + .hk)s + j2(.hh - l)p 2 
+ 0(p~s,s 2 ), (2.13) 

where qo, hq, tuq, vq, n 3 , etc. are the quadratic mixing 
coefficients entering the full nonlinear scaling fields: see 
I(1.4)-(1.6); in addition, one discovers that the finite-size 



L~ d correction terms in the scaling fields — see (1.8) — 
enter only with the quadratic mixing coefficients. Like- 
wise we obtain 

s = k + (v + k n )fi + (n 3 + 2k m Q )p + (2r + k n 3 )t 
- (fci +.nko)p+ (1 - ji/co)S + OG5 2 ,/5S,S 2 ), (2.14) 

where, again, we have retained only the leading terms 
needed later: further terms are given in K(4.29)-(4.30). 

Similarly, the generalized susceptibilities defined in 
1(2.28) are useful here: one finds 

Xhh = {d 2 P/dh\ = U 2 L L^(d 2 y Y)/p c , (2.15) 

and likewise for \ht an d Xtt- The basic number fluctu- 
ation or reduced susceptibility xnn — {d 2 p/dp 2 ) t can 
then — see 1(2.29) and K(4.33) and Appendix F — be 
expressed as 

PcXnn = e 2 U 2 L^(d 2 y Y) 

-3j 2 e 1 UlL^-^(d 2 y Y)(d y Y)/p c 

- 2e 1 e 3 U L D L L^-^/"{d x d y Y) + ■■■, (2.16) 

where only the leading terms have been displayed while 
the constants are 

ei = 1 - .72, e 3 = h + ii, (l = 1); (2.17) 

sec 1(2.30) and 1(3.22). This result is needed to study 
the fc-loci in finite systems: see Sec. III. A. The Q-loci, 
taken up in Sec. IV. A, demand the higher-order analogs. 

D. Chemical potential at (T c , p c ) 

Before turning to the various critical loci and their 
finite-size behavior, we address a rather special question 
which turns out to be interesting since its answer, as 
mentioned in the Introduction, opens an opportunity to 
determine via precise simulations, the presence or ab- 
sence of finite-size dependence in the scaling fields. In 
a finite grand canonical ensemble at temperature T the 
chemical potential p must be adjusted to achieve a spec- 
ified density: but the resulting value will depend on L. 
Accordingly we ask: "How does the finite-size chemical 
potential, say p c L = pl(T c ,Pc), needed to achieve the 
bulk critical density, p c , at the critical temperature, T c , 
approach p^ = p c T 

To attack the problem we first determine the scaling 
fields at T = T c and p = p c , i.e., t = and p = p c = 1. 
Recalling that l = 1 [1(3.22)], the relation (2.13) for the 
density p then yields 

= (1 - h)f> - (h + 3i)s - j 2 (l - j 2 )p 2 + • • • , (2.18) 

where we have neglected the "background" terms in p 
and p (arising from the quadratic mixing coefficients) 
and may check later that they yield only higher-order 
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corrections. (Note that § ~ L _ ( 1_Q )/ l/ dominates L~ d 
since d > 2 > (1 — a) /v.) By appealing to (2.11) and 
the scaling function expansions (2.9) and (2.8) this can 
be re-expressed as 

2(l-j2)U L {Y ° 2 + Y^y c L4 + ---}y L 

-(h + h)D L LV- 1+a V v [Y? + -..]« 0. (2.19) 

From the definitions (2.3) of j/l and yhk we thus find that 
when p = p c at t = the ordering field obeys 

h « ap/lP-"^" = a^/L^- 1 ^, (2.20) 

where the omitted correction factor includes L~ 8i l v and 
j^-i/u+d-d as i eac [i n g contributions, while 

a„ = (h + ji)D L Y? /2{l - 3 2)U 2 A (2.21) 

Note also that even in the absence of pressure mixing 
(i.e., ji = j 2 = 0) the contribution of p to t, via l x ^ 0, 
ensures that h does not vanish (as it would identically 
in a symmetric system); instead h decays with a leading 
exponent d + (7 — l)/v of value about 3.38 for the d = 3 
Ising universality class. 

Finally, at t = the relation (1.4) for p with the added 
term —s 2 /L d , and (1.8), leads, in linear order, to 



p{dXNN/dfi) T = Kxnn) 2 



(3.1) 



p, = h + j 2 p + s 2 /L d 
= h + j 2 (p + A + s /L d ) + s 2 /L d . 



(2.22) 



On using (2.2) for p at xl ~ J/z, ~ this may be solved 
to yield 

Al = K^cPcSi) - Aic]/fc B r c , 

= a L /L* + a p /L d + „/(l - j^+^Z" + • • • , 

(2.23) 

where the new amplitudes are 

= (S 2 + J2S0)/ (1 - .72), ftp = j2*oo/Pc(l - .72)- 

(2.24) 

Evidently, if d < d and J2S0 and S2 do not both vanish, 
the dominant behavior arises from the L-dependence of 
the scaling fields. If pressure mixing is absent (or negli- 
gible) the last, most rapidly decaying term in (2.23) will 
be controlling. 



III. MODIFIED-SUSCEPTIBILITY LOCI IN 
FINITE SYSTEMS 

A. Asymptotic expressions 

The /c-modified-susceptibility loci or, for brevity, the k- 
loci are defined by the isothermal maxima of \^ = XI 1 P k 
and so satisfy 1(4.32), namely, 



We aim to solve this equation asymptotically near criti- 
cality, first, to obtain fS k \t;L), i.e., the finite-size fc-loci 
in the (p,, T) plane, then p( k \t;L), and, finally, p {k) {t;L), 
the locus in the (p, T) plane which is of most practical 
interest. 

The required third order susceptibility, \n 3 = 
(dxNN/dp)T, can be obtained by differentiating (2.16) 
with respect to p at fixed t. This entails the derivatives 

(dx L /dp) T = D L I}I V {-1 X - hp +•••), (3.2) 

(dy L /dp)T = U L L A /»(1 - j 2 p + •••), (3.3) 

which follow from (2.3), (1.3), (1.4), and (2.12). On using 
(2.13) for p this leads to 

PcXN s = e\UlL^y»{dlY) 3 2 e\p?UiL*</» 

x [A(d 3 y Y)(d v Y) + 3(d 2 Y) 2 } 

- ie\e 3 UlD L L^l v {d x dlY) + ■■■, (3.4) 

where we recall (2.17) for e\ and e%. Using the expansions 
(2.6) and then (2.8), for the scaling functions Y K (xL,yL), 
yields, after some algebra, the defining equation (3.1) in 
the form 



24 ei y ° 4 + 24e 1 y i > L + 2Ae 1 Ul 4 Y^ ) L- B l" 

- (3.72 + ke x )e x p- 1 \J L L- si l v [2Y ° 2 + 2Y 1 ° 2 x L 

- Ze z {D L /U L )L^-^^ [2Y° + 2Y 2 \x L + • • 
+ ••• = 0. 



VL 



(3.5) 



With the aid of (2.3) the scaling field h can hence be 
written in terms of L and t as 

h = ^(3j 2 + ke 1 )/p c Y 4 L ( - 2 - a y^2Y° 2 + 2Y* 2 D L tL 1 ' v 

+ 2Ui 4 Y$L-^ + • • -] 2 - Ul 4 Y fh/Y ° 4 L^ 

+■■■• (3.6) 

In order to solve this equation for p as a function of 
L and t, we first write p in terms of fi, t, and L by us- 
ing the finite-size scaling equation (2.2). The expansions 
(2.8) for Y(xl, ■ ■ •) can then be employed and on solving 
iteratively for p we obtain 

PcP = Pc(k t + fi + s L- d +•••)+ Y^L- {2 -^' v 
+ D L Y? Q [(l-j 1 k Q )t-(h + rirtL-^l" 
+ C/£ 4 r ( 4) L-( 2 -«+^ + .... (3.7) 

Rewriting (3.6) yields the reduced chemical potential, p, 
in a similar form from which p may be eliminated using 
(3.7). Solving for p iteratively as a function of t and L, 
finally yields the finite-size /c-loci in the (/x, T) plane as 

p^(t;L) = [p^(T;L)-p c ]/k B T c 

= p[ k) t + (S2 + hs )L- d + M[ k) L-V- a V v 

+ Mi k) L-( 2 - a+ W» + M^iL-^/" 



+ 



(3.8) 
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where the Mf> vary linearly with fc and are given explic- 
itly in K(4.53)-(4.54) while /l ik) = (fei +j 2 k )/(l -j 2 ), is 
actually independent of k and equal to /t CT; i which was 
defined in 1(3.16) as the (reduced) slope of the phase 
boundary H a (T) at T = T c . Notice that, owing to the 
hyperscaling relation, the L~ d term here dominates the 
universal scaling contribution, _L _ ( 2_Q )/ l/ = L~ d , when 
d < d. 

Substituting (3.8) in (3.7) yields the fc-loci in the (p, T) 
plane as 

p^(t;L) = \p^(T;L)-p c }/p c k B T c 

= p[ k) t+[(l+h)s + s 2 }L- d 
+ (M[ k) +Y^)L- {2 -^' v 

+ (Mf ) + c/£ 4 r ( 4) )L-( 2 -«+^ 

+ (M 3 W + D L Y° r)tL-^/ v + • • • , (3.9) 

where p[ k ^ = fco + fi^ is also independent of k and equal 
to p CTi i [ see 1(3.12)] while 

t = 1 - jiko - (h + ii)(fci + j 2 fc )/(l - .72), (3.10) 

which, in fact, has the same value as r in 1(3.14). 

To obtain the fc-loci in the (/?, T) plane, we now sub- 
stitute (3.8) and (3.9) into the scaling fields h and t to 
find 

VL = W32 + k ei )U L (Y^f/Y ° 4 L^ 

x[i + 2c/£ 4 y ( 2 4) /r °2^ /!/ 

+ 2D L Y 1 ° 2 rtL 1 ^/Y ° 2 + •••], (3.11) 
x L = D L rtL 1/v + • • • , (3.12) 

and thence can express the generalized densities, p and 
s, in (2.11) in terms of L and t. Finally, from (2.13) we 
obtain the desired fc-loci in the (p, T) plane as 

p {k \T-L)/ Pc = 1 + B (k) L- 2t) /» + C {k) L-^l v 

+ B (k) L -W+o)/" + ... + A ( f ) t + --- 

+ A 2 k) L- d + B {k) L-« 3+e ^l v + ■ • • , (3.13) 

where the leading coefficients are 

B (k) = (1 - j 2 )(3j 2 + k ei )U 2 L (Y 2 f/S P 2 X 4 , (3.14) 
C[ k) = -(h+j 1 )D L Y 1 ° /p c , 

Bi k) =3B[ k) Ul 4 Y ^/Y ° 2 , (3.15) 

A ik) =v Q + n 3 + (2q + n )fi {k) + (n + 2m )p[ k) , (3.16) 

while A 2 k ^ and , which also entail the nonlinear 
scaling- variable coefficients v , n , m , q , ■ ■ ■ [see 1(1.4)- 
(1.6)], are given in K(4.62). 

Note that the coefficient A\ of the leading analytic, 
L-independent term actually coincides with A\ in 1(3.26) 



which is the amplitude of the linear t term in the co- 
existence curve diameter. Furthermore, the contribu- 
tion from the finite-size corrections to the scaling fields, 
i.e., the L~ 3 term in (3.13) is dominated by L~ 2/3/V , 
L -(i- a )/u and L -(2p+e)/u tcrms ( provided d > 2 ). When 

(k) 

T = T c , the analytic, ^-independent part of p c vanishes. 
The leading correction then decays as i -2 ' 3 /^ with an 
amplitude that varies linearly with k; this is followed by 
an L^ 1- "'/ 1 ' term whose amplitude does not depend on 
k. As mentioned in the Introduction, the leading ampli- 

(k) 

tude, B\ , vanishes, in fact, when k assumes the "opti- 
mal value" fc pt = — 3j2/ei = 3i? M , where the Yang- Yang 
ratio is defined in Ref. [8] and I Sec. III.E. This value 
coincides with the one obtained in 1(4.37) for the ther- 
modynamic limit when it should describe the particular 
fc-locus that approaches the critical point "most directly" 
in the (p, T) plane. 

B. Finite-size fc-loci: behavior and applications 

The near-critical behavior of the finite-size fc-loci for 
the hard-core square-well (HCSW) fluid and for the re- 
stricted primitive model (RPM) electrolyte is illustrated 
in Figs. 1(a) and (b), respectively. The results shown are 
based on simulations in periodic cubical boxes of dimen- 
sions L* (= L/a, where a is the hard-core diameter) upto 
13.5 and 12, respectively [7,19]. The limiting (L — ► 00) 
behavior for the same models is shown in Figs. 1 and 2 
of I (while results for a van der Waals fluid are shown in 
I Fig. 3). The differences between the HCSW and RPM 
are quite striking: for the former a value of fc op t close to 
zero or even somewhat negative is suggested while for the 
RPM one might conclude fc opt ~ 0.8. These (inevitably 
rather uncertain) estimates correspond surprizingly well 
via fcopt = 3-R^ with more recent (quite independent) 
estimates for the Yang- Yang ratio of —0.044(3) and 
+0.26(4) for the two models [25]. 

The result (3.13) shows that the density estimated at 
T = T c on the fc-locus, namely, p^(T c ;L), approaches 
the bulk critical density, p c , as, in leading order, L~^, 
with if) = 2/3 jv provided the pressure mixing coefficient 
j'2 does not vanish. For [d = 3) Ising-type criticality this 
predicts ip ~ 1.03 whereas for a classical system tp = 2. 
If ]2 (and, hence, vanishes or is numerically small, 
the next leading term in (3.13), varying as L _ ( 1_Q )/ Iy j 
becomes dominant. The exponent ip = (1 — a)/v then 
takes the value 2 for classical criticality but ~ 1.41 for 
(d = 3) Ising systems. 

If a reliable estimate for T c is known — we indicate 
below [in Sec. IV. C] how this may be found by using the 
Q-loci — these results can be used in simulations to ob- 
tain convincing, unbiased estimates of the critical density 
p c . By "unbiased" we mean that prior knowledge of the 
critical universality class is not required. One effective 
strategy is implemented in Fig. 2 for the HCSW fluid 
where p*{L*) = pW(T = T c ;L*)a 3 has been plotted for 
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FIG. 1. The fc-loci in the (p,T) plane for (a) the hard-core 
square- well fluid with, from the right, k — 0, 0.25 and 1 where 
the system sizes L* used in the figure are 5, 6, 7.5, 9, 10.5, 12, 
and 13.5 (measured in units of the hard-core diameter, a = a) 
[7]; and (b) the restricted primitive model electrolyte with 
k — 0, 0.5 and 1 where the system sizes shown are L* = 6, 7, 
8, 9, 10, and 12 [19]. Note that p* = pa 3 while the reduced 
temperatures T* are defined in Refs [7] and [19] and in Sec. 
V below. 

k = and k = 1 vs. 1/ L*^ for trial values of the exponent 
ip varying from 1 to 2 (which range encompasses both 
the classical and (d = 3) Ising universality classes). For 
these plots the HCSW estimate T* = k B T c /e ~ 1.2186, 
obtained in Sec. IV. C below, has been used. It turns out, 
however, that p^(T;L*) is rather insensitive to T ~ T c 
so that essentially the same results are obtained if the 
original Orkoulas et al. [7] estimate (which is about 0.06% 
lower) is used instead. (Note that this insensitivity is not 
realized in the RPM !) 

The straightest plot for k = 1 [in Fig. 2(b)] corresponds 
to -0 — 1-0 which is consistent with Ising behavior (as 
expected). However, the k = plots in Fig. 2(a) are 
straightest for tp = 1.4-1.7: this is also consistent with 
Ising behavior provided (as seems to be the case) the 
value of j2 is small. Together these plots suggest a critical 



FIG. 2. The scaling behavior of p* c {L*) at T = T c for (a) 
the k = locus (solid lines and squares) and (b) the k = 1 
locus for a hard-core square-well fluid [7] for trial values of 
the exponent ip: (i) 1.0, (ii) 1.2, (iii) 1.4, (iv) 1.7 and (v) 
2.0. The dotted lines and open squares in part (a) derive from 
the Q-loci: see Sec. IV. A. 

value of p* in the range 0.3065 to 0.3080. To improve the 
possibilities for extrapolation, the k = data are com- 
bined with data for k = 0.25 and 0.1 in Fig. 3 and plotted 
vs. A/(L* + l*)^ , where A is merely a convenient scale 
factor while the "shift" I* has been introduced to allow 
(approximately) for the anticipated higher order correc- 
tions. From this figure, we estimate p c for the HCSW 
fluid (with interaction range b = 1.5a [7]) as 

p* = p c a 3 = 0.3068 ± 0.0007. (3.17) 

This value agrees well with Orkoulas et al. [7] who found 
p* = 0.3067 ± 0.0004. By the same approach Luijten et 
al. [19] estimated the critical density of the RPM elec- 
trolyte but only to the rather lower precision of ±3% 
which, however, should be more reliable than other, less 
systematic and biased methods. 
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FIG. 3. Estimation of the critical density for the HCSW 
fluid by extrapolation to L — > oo. The upper solid sym- 
bols derive from the k~0 locus with, from the right, 
(1>X,A) = (1.0,-1.5,0.7), (1.2,-0.5,1.0), (1.4,0,1.4), 
(1.7,1.0,2.4), (2.0,1.5,3.5). The central crosses, from 
the fc = 0.1 locus, have (tp,l*,A) = (1,0,1). The 
lower, open symbols are plotted with, from the right, 
(ip,l*,A) = (1.0,0.5,1.0), (1.2,2.0,1.7), (1.4,3.0,2.7), 
(1.7,4.5,5.6), and (2.0,6.0, 10.5). 



IV. BEHAVIOR OF THE Q PARAMETER AND 
Q-LOCI 

Some time ago Binder [20] introduced the dimen- 
sionless, finite-system moment ratio, Ql{T;(p)l) = 
(to 2 )|/(to 4 )l, defined in a grand canonical ensemble 
with m = p — (p)l, and showed how, in simulations of 
symmetric systems (where p = p c is known), it was par- 
ticularly useful in locating the critical temperature pre- 
cisely. Specifically, plots of Ql(T;p c ), evaluated on the 
(symmetric) critical isochore at values of L increased in 
steps by increments Ai, display successive intersections 
at temperatures, say, Tq L (L), that rapidly approach the 
limiting, critical temperature T — T c . At the same time 
the intersections define a unique and universal critical 
value [21-23] Q c = hin^^oc Ql(T c ; p c ). However, the 
obvious difficulty in attempting to adapt this approach 
to a nonsymmetric fluid system is that the critical den- 
sity is not known; nor, in fact, even if p c were known, 
is it clear that the critical isochore would be the most 
appropriate locus on which to examine the temperature 
dependence of Ql- Indeed, we will see from our study of 
Ql(T; (p)l) for general systems that is presented here, 
that the locus p = p c , even if known, would not normally 
be optimal! 

To make progress as explained in the Introduction, 
we define, following [19], the Q-loci, pq(T;L), via the 
isothermal maxima of Ql(T; (p)l) where, it is worth 



reemphasizing, (-)l denotes a grand canonical finite-size 
average in which p is chosen to yield the desired values of 
the mean density {p) l (which, of course, is distinct from 
what might be considered for a canonical system in which 
p = N/V is directly controlled and does not fluctuate). 
As seen in Fig. 4, for the HCSW fluid and the RPM, 
the ratio Ql at fixed T displays a unique maximum vs. 




FIG. 4. The moment-ratio parameter Ql(T;p) vs. p at 
fixed temperatures (a) for the hard-core square-well fluid at 
L* = 10.5 (from the top, T* = 1.0, 1.1, 1.15, 1.2, 1.22, 1.25, 
1.3, 1.35, and 1.4): note that T* ~ 1.2179 [7], and (b) for 
the restricted primitive model electrolyte at L* = 10 (from 
the bottom, the solid-lines are for 1/T* = 13, 15-19, 19.5, 20, 
20.5, and 21); the dashed line is at T* ~ 0.050 [19]. 

density so that pq(T;L) is well defined. In more com- 
plex models with, e.g., more than one critical point, the 
loci will presumably display separate branches or more 
complex topology; but our concern here is with the be- 
havior of the loci near criticality as L — > oo, first in the 
one-phase region above T c , then through the two-phase 
region below T c . In the following section we illustrate the 
explicit use of these results in simulations. 
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A. Q-loci above criticality 

As observed originally by Binder, thermodynamic den- 
sity fluctuations in a single-phase region of the phase 
plane should follow a Gaussian distribution when L — > oo 
so that Ql(T; {p)l) for T > T c should tend to the con- 
stant value 4 as L increases. In practice, as illustrated 
in Fig. 5, the approach at fixed T is nonmonotonic and 
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FIG. 5. Variation of the moment ratio, Ql{T;p), 
with increasing size for a hard-core square-well fluid 
at T* = 1.300 ~ 1.0674 T c * [7]; the system dimensions are 
L* — 5, 6, 7.5, 9, 10.5, and 12. The horizontal solid line rep- 
resents the single-phase limit Qoo = §■ 

entails a progression of the Q-locus to an apparently well 
defined limit p%(T). 

To estimate the asymptotic behavior of pq(T;L) we 
may follow the strategy used in studying the /c-loci. 
First, in terms of the generalized susceptibilities XN k = 
(d k p/ dp k )r with p = p/k^T and p = p/ksT, note that 
Ql is equivalent to V(xnn) 2 /xjv 4 - Thence we find 

\ dp J T {xn±V 

from which, since (p) l increases monotonically with p at 
fixed T, one sees that the pq locus satisfies the equation 

2Xn*Xn* - XawXw 5 = 0. (4.2) 

Here we have employed the reduced susceptibilities XNNi 
Xn 3 , etc., introduced in (2.16) and (3.4). From (3.4) we 
then obtain 

PcXm = e\UlL^/» [{d*Y){dy L /dp) T 

+ (d x d 3 y Y)(dx L /dp) T ] 

- help-'UiL 2 ^ [A(d 4 y Y)(d y Y) 



+ io(d 3 y Y)(d 2 Y)] (d yL /dp) T 

j 2 e 3 1 p- 1 U 4 L L 2 ^ [A(d x d 3 y Y)(d y Y) 

+ A(d 3 y Y)(d x d y Y) + 6(d x d 2 y Y)(d 2 Y)] (dx L /dp) T 

- 3e 2 e 3 U 2 L D L L^ +i y" [(d x d 3 y Y)(dy L /dp) T 

+ (d 2 d 2 y Y)(dx L /dp) T ]+---. (4.3) 

Using (3.3) for (dyL/dp)T and, in that result, (2.13) for 
p yields 

P.XN* = etuiL^y^Y) Sj 2 e\p-'UIL^+^ 
x[(d*Y)(d y Y) + 2(d 3 Y)(d 2 Y)] 

- ^e\e 3 U 3 L D L L^ + ^' v {d x d 3 y Y) + ■■■. (4.4) 
Similarly, after some algebra, we obtain 

PcXivs = e\UiL^)'»{dlY) - j 2 e\p-'UlL 2 ^/» 
x[6(d 5 y Y)(d y Y) + 15(d 4 y Y)(d 2 y Y) + lQ(d 3 Y) 2 } 

- 5eie 3 UiD L L^+ 2A +^^(d x d^Y) + ■■■. (4.5) 

Now we may use the leading approximation (2.9) for 
the scaling function Y and substitute the expressions 
(2.16), (3.4), (4.4) and (4.5) into the Q-locus equation 
(4.2). This then reduces to 

[40O 2 - 5Y ° 2 Y l]y L + 3j 2 (Y° 2 ) 2 Y Q ° 4 U L /p c L^ + ■ ■ ■ = 0, 

(4.6) 

where, for brevity we have displayed only the leading 
terms; this, in turn, is readily solved to yield on the 
Q-locus as to 

(i ^_.nYQ 3(y " 2 ) 2 y 4 t/ L / Pc 

To obtain the density, p, we appeal to (2.13) and use 
(2.11) for p and s; then, with (2.6) and (2.8) for the 
scaling functions, and using (4.7) for y^, we finally obtain 
the Q - l° cus explicitly as 

PQ (T; L)/p c = 1 + BqL- 2 ?/" + C Q L-^I» 

+ A Q t+---, (4.8) 

where the leading coefficients are 

B Q = -2j 2 (l-j 2 )Y Q Y° 2 U L /p c , 

C Q = -(h+j 1 )Y 1 ° D L /p c , (4.9) 

while Aq is equal to = A\, the (reduced) slope of 
the coexistence curve diameter as given in (3.16) and 
1(3.26). Note that the leading amplitude, Bq, vanishes 
when j2 = 0. As an explicit example, we present the 
Q-locus for a hard-core square- well fluid [7] in Fig. 6. 
Evidently, the loci both above and below T c approach 
the critical point when L — > oo. 

As illustrated by the open squares in Fig. 2(a) above, 
the evaluation of pq(T;L) at T = T c can be used to 
provide unbiased estimators for the critical density, p c , 
that, in fact, resemble quite closely the sequence provided 
by the k = loci: see also Fig. 5 in Ref. [19]. 
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P* 

FIG. 6. The Q-loci in the (p, T) plane for a hard-core 
square-well fluid. From the right, the simulation box dimen- 
sions are L* = 5, 6, 7.5, 9, 10.5, 12, and 13.5. The estimated 
critical point shown is {pt,T*) = (0.3067, 1.2179) [7]; the solid 
dots represent the estimated coexistence curve diameter. 

B. Modified or Q (fc) -loci 

It is instructive to define a modified Q parameter, as 
for in Sec. III. A, via 

Q^(T;{p) L ) = Q L {T;{p) L )/{p) k L . (4.10) 

The modified or Q^-loci are then defined by the points 
of isothermal maxima of in the (p, T) plane. In 
terms of the reduced susceptibilities, the equation for the 
locus Pq\T;L), becomes 

p[2xn 3 Xn* ~ XnnXn 5 } = kx 2 NN XN*, (4.11) 

which extends (4.2). The extension of (4.6) gains the 
term \k{\ - j 2 ){Y§ 2 ) 2 Y§ i U L L-l 3 / v / p c on the right hand 

side. Finally, p^ (T; L) is represented by the same ex- 
pression (4.8) [for k = 0] except that Bq must be replaced 

b y b q ] s iven b y 

j 2 B™ =B Q [j 2 -\k{l-j 2 )] . (4.12) 

This coefficient vanishes when k = kQ = %j 2 / (1 — j 2 ) = 
— 9Rfi which may be contrasted with the "optimal" k- 
locus specified by k opt = 3i? M (see Sec. III. A). Note that 
the coefficients Cq and Aq in (4.8) do not gain any k de- 
pendence although various higher order coefficients will, 
in fact, depend nonlinearly on k. 

C. Behavior of Q in the two-phase region 

At fixed T < T c the phase transition in the thermo- 
dynamic limit is of first-order character with a jump 



in density from p-(T) to p+(T) as p increases through 
the phase boundary, p a (T). Finite-size scaling theory 
has been extended to first-order transitions [20,37-40] 
although the main focus previously has been on the de- 
pendence as a function of the field h cx p — p a (T). Here, 
motivated by the requirements of simulations, we will en- 
quire more closely into the variation with the density, p. 
From this perspective, the crucial feature is that when 
p ~ p a the grand canonical equilibrium distribution 
function, Pl(p', Pi T), exhibits two peaks located at densi- 
ties near p_(= p vap ) and p+(= pn q ). For sufficiently large 
L these peaks can be represented as Gaussians [20,39,40]. 
Inside the two-phase region one may also need to con- 
sider the surface free energy associated with interfaces 
that separate domains of coexisting phases [24,37,41,42]. 
However, for regularly shaped domains (such as periodic 
cubes or fixed-shape parallelepipeds) these contributions 
enter only as exponentially smaller corrections, so they 
are not considered here. In the case of general fluids the 
density distribution Pl{p) has no symmetry: thus for 
large L in a d-dimensional system of fixed regular shape 
with periodic boundary conditions, we will accept the 
form [24] 

P L (p: p, T) « C L { X Z 1/2 exp[-/3(p - p_) 2 L d /2 X -] 

+ X+ 1/2 cxp[-[3(p-p + ) 2 L d /2 X+ }} 
xexp[/3p(p- p a )L d ], (4.13) 

where (3 = 1/ksT, while Cl(p, T) is a normalization con- 
stant, and the x±(T) are the infinite- volume susceptibil- 
ities [defined via x = (9p/dp)T] at p = p±(T)±. This 
distribution has been set up so that when p = p a both 
Gaussians contribute to Pl (p) with equal weight [43] . 

To simplify subsequent expressions let us introduce the 
basic, dimensionless ordering field 

h=[p-p a (T)}/k B T, (4.14) 

and the average and difference densities and susceptibil- 
ities 

P(T) = \( P+ + P-) and p (T) = i(p+-p_), (4.15) 
X(T) = i( X+ +X_) and Xo (T) = ±(x+ - X-)- (4.16) 

Note that xo vanishes identically in a symmetric sys- 
tem. For further convenience here we also define the 
augmented field-dependent densities 

P + = P + X\h, Pq = Po + Xoh, and p ( Q h) = p + \xo h - 

(4.17) 

By replacing the summation over discrete density val- 
ues, p = N/V > 0, by integration over p and extending 
the lower limit to p = — 00 (which will entail only an 
exponentially small error for large L), we may compute 
(p)l and the moments This yields 
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(p) L (p,T) w p+ + P +t & nHhp i i l) L d ). (4.18) 

Note that when h = 0, or p = fi a (T), we have (p)l ~ 
p(T), i.e., the coexistence curve diameter. Likewise we 
find 



(m 2 ) L (p,T)^f + f 1 /PL d , 

(to 4 )^, T) « / 2 + / 3 //?L d + ,A//3 2 i 2d , 



where, with 



A/? =(p)l~ P + and T = tanh(/ip[ ) ' i) J L d ), 



(4.19) 
(4.20) 



(4.21) 



From these results it is evident that Ql((p)) is a ratio of 
two polynomials of fourth order in (p) but quadratic in 
L- d . 

To examine the two-phase behavior of Ql in the ther- 
modynamic limit, let us define the scaled deviation from 
the coexistence diameter, p(T), via 



y = (p- p)/po, 



(4.26) 



so that y = ±1 for p = p±(T). In the first instance we 
may then, as in [19], set p = p a (or h = 0) before allowing 
L — > oo. As observed after (4.18) we then have (p)l — > 
(p)oo = P and T = in (4.21)-(4.25). If nonetheless, we 
identify (p)l in (4.21) as p in (4.26) and evaluate (m 2 )^ 
and (m 4 )oc accordingly one is led to 



g^(r; /5 ) = l-4y 2 /(l + 6 y 2 + 2/ 4 ), 



(4.27) 



which, apart from the superscript a which indicates the 
limiting procedure adopted, is the result quoted, mislead- 
ingly, in [19]! Indeed, this can only be the correct limit 
of Ql(T; (p)) when T < T c if y — 0, i.e., on the diameter. 

To obtain the true limiting behavior for — 1 < y < 1, 
one must first notice that for (p) l to approach a general 
value in the interval the thermodynamic limit 

must be taken with hL d in (4.18) approaching a finite 
value that yields (p)z — > p for the desired value of y. 
This corresponds, in fact, to T w tanh(/ip ^ d ) ~ 2/ an d 
then yields — see also [24] — the limiting moments 



(m 2 )oc = Po(l - y 2 ), 
(m 4 ) oc =p 4 (l-y 2 )(l + 3j/ 2 ), 



(4.28) 
(4.29) 



both of which, perhaps surprisingly, vanish linearly on 
the phase boundary, i.e., as y 2 — > 1 — . Equally, then [24] 

Q 00 (T;(p)) = (l-y 2 )/(l + 3y 2 ) (T < T c ), (4.30) 



vanishes linearly on the phase boundary. On the other 
hand, Qoo{T < T c ) takes its maximal value, namely 1, 
on the coexistence diameter (y — 0). Indeed, the corre- 
sponding approach of the Pq{T; (p)) loci below T c to the 
diameter is evident in Fig. 6. 

To give a graphic impression of the limiting behavior of 
Ql(T; (p)) we display in Fig. 7, plots constructed using 
(4.30) and the coexistence curve data for the hard-core 
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FIG. 7. The behavior of the limiting moment ratio, 
Qoc{T;p), vs. p at fixed temperatures below T c for the 
hard-core square- well fluid [7]. The solid line, dashed 
line, and long-dashed line are for T/T c ~ 0.82, 0.90, and 
0.985, respectively. The cross is at the critical point 
(r*,p*)~ (1.218,0.306) [7]. 

square- well fluid [7] at various temperatures below T c . 
In addition we have indicated by a cross the anticipated 
Ising critical point value, Q c — 0.6236(2) [22,23], that we 
also verify independently below. The horizontal line at 
Q = 3 describes the limiting single-phase value. 



D. Scaling of Ql({p)) near coexistence 

In the one-phase region outside the coexistence curve, 
i.e., for y 2 > 1 [sec (4.26)] the result Qoo — | should be 
recaptured by the analysis based on (4.13); indeed, the 
results (4.18)-(4.25) do confirm this. Thus for h nonzero 
and L — > oo, the expression (4.18) yields 

(p) L ~P±+ X±h - 2(po + X vh)e- 2hpaL \ (4.31) 

where the + or — corresponds to h ^ 0. On substitution 
in (4.22)-(4.25) the L-independent terms in (m 2 )^ and 
(to 4 ) l cancel identically leaving 

(m 2 ) L - X±/PL d + 0(e- 2h ?° L \ (4.32) 
and, similarly, (to 4 )l 3(to 2 )| / , yielding finally 
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Ql(T) 



l + 0(e~ 



2hp L c 



(4.33) 



for T < T c and h nonvanishing (but not too large). 

Evidently, in the thermodynamic limit, Qoo(T; (p)) 
vanishes as p approaches p + or p_ from the two-phase 
region and then jumps discontinuously to ^ on entering 
the single-phase domain. This behavior as L — > oo can be 
seen clearly in grand canonical simulations as illustrated 
in Fig. 8 for the hard-core square- well fluid [7]. The pre- 




FIG. 8. Behavior of Ql(T;p) for a hard-core square-well 
fluid at T/T c ~ 0.944 [7]. The thin lines represent simulation 
data for L* = 5, 6, 7.5, 9, f0.5, and 12, while the thick line is 
the prediction for L — oo [scaled to the estimated values of 
pX(T) and p*_(T)}. 

dieted limiting behavior is approached rather rapidly at 
the selected temperature, namely, ~ 5% below criticality. 
However, closer to T c and for the RPM the convergence 
is much slower and less regular as seen in Fig. 9 which 
reports simulations ~1.5% below the (estimated) criti- 
cal points. In all cases — as follows from previous the- 
oretical and simulation-based observations [24,40,41] - 
the plots of Ql(T; (p)) display rounded, but increasingly 
deep and sharp minima outside, but approaching, the 
coexistence curve as L increases. However, the strongly 
asymmetric and relatively slow approach of the RPM to 
the limiting behavior is striking. Nevertheless, it turns 
out that by tracking these minima and suitably extrapo- 
lating them on the basis of the present theoretical founda- 
tions, remarkably precise estimates of the density jump, 
2po (T) = p+(T)-p_(T), and of the diameter, p(T), can 
be obtained for both models [25]. 

In order to understand the minima better let us, for 
simplicity, consider the symmetric case where x+ = X- 
so xo = in (4.13)-(4.16). After some algebra we obtain 
from (4.19)-(4.25) the expression 
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FIG. 9. Simulation data for Ql(T;p) (a) for the hard-core 
square-well fluid at T/T c ~ 0.985 using the same box sizes L* 
as in Fig. 8; (b) for the RPM electrolyte at T/T c ~ 0.986 for 
L* = 5-10 [19]. The thick lines represent predictions for the 
limit L — oo. 

n (T ] = [X + (l-T 2 )} 2 

vm -.p) 3A-2 + 6A-(i-r 2 ) + i + 2r 2 -3r 4 ' 

(4.34) 



where T(hp L d ) was defined in (4.21) while 
X(T,h;L)=x(T)/p 2 Q (T)k B TL d . 



(4.35) 



When L oo so X and h -» with T 2 y 2 < 1, 
the previous result (4.30) is recaptured; on the other 
hand, when L — > oo with h fixed and nonzero, one has 
T 2 -> 1 and (4.33) is matched. A plot of Q L vs. y gen- 
erated from (4.34) [see K Figs. 4.9 and 4.8] quite closely 
mirrors, except for its precise y <4> — y symmetry, the sim- 
ulations for the HCSW fluid shown in Figs. 8 and 9(a): 
indeed, the HCSW fluid does not deviate drastically from 
overall symmetry even though it displays some pressure 
and chemical potential mixing (as discussed above). 
For finite L, (4.34) predicts two minima that satisfy 



T± =± (l + 2X) 1/2 /{l + 3X) 1/2 , (4.36) 
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Qmin(T' : L) 



X{2 + ZX 2 ) 



2\2 



4 + 18 A" + 36^ 2 + ' 
x/p 2 k B TL d + 0(e- 2h r° Ld 



(4.37) 



Thus Q m in(L) approaches zero, the limiting value at co- 
existence, as L~ d . On the other hand the positions of the 
minima approach p + and p_ when L — > oo. In order to 
find the corresponding L dependence, we first determine 
h± from (4.36) and (4.21) obtaining a (\nL)/L d varia- 
tion. From (4.18) we thence find the density minima at 

ftinW L) = p±(T) ± 2p (T)B Q (T)L- d 

x[\n(L d /B Q )-l + 0(L-% (4.38) 



where the scaling amplitude is 

B Q (T) = k B Tx{T)/Apl(T). 



(4.39) 



Since this result has been derived only for the symmetric 
case (although it has wider validity [25]) we may replace 
X by x+ = X-\ ^ is also useful to recall that 2po = 
P+-P- = Apoo(T) [sec (1.11)]. 

Our discussion of Ql(T; (p)) below T c has, up to this 
point, been confined to fixed T and, then, to large enough 
L. On the other hand, when t = (T — T c )/T c -> 0- the 
basic thermodynamic properties entering the expressions 
for Ql(T; (p)) and for the minima and their locations 
will display their standard critical behavior, specifically, 
Po ~ \t\ fi , X ~ \t\-\ while Xo ~ \t\^ [sec 1(3.41,3.42)]. 
Beyond that, however, the divergence of the correlation 
length, namely, £ ~ a/|t|", implies that each variable L 
appearing in the formulas above should, when t — > 0— , be 
associated with a factor \t\" . However, the analysis based 
on the two-Gaussian form (4.13) implicitly assumed that 
w = L/£ ~ L*\t\" was large [24,37-42]: thus when t -> 
0— , we may not simply substitute the expected powers 
of t in to the expressions so-far derived. On the other 
hand, the full scaling expression for Ql implied by the 
basic scaling ansatz (2.2), namely, 



Ql(T;p) w Q(x L ,Y L ,y L4 , y L5 , ■ ■ •), 



(4.40) 



must reproduce the expressions obtained here when 
w = L/£ ~ \xl\ v — > oo [sec (2.3)]. This means 
that although we cannot hope to derive theoretically 
an explicit general expression for the scaling function 
Q(x,y, •••), or even the scaling forms for the reduced 
minima, p^ in {T; L)/p (T), we have in essence obtained 
exact information about the corresponding scaling behav- 
ior! It thus transpire, as shown in [25], that by starting 
at a temperature below T c where £(T)/a = 0(1), simula- 
tion data at increasing T can be used to generate the ap- 
propriate scaling functions for p^ in (T; L) and p^ n {T; L) 
and thereby also obtain precise estimates for po{T) and 
p{T), i.e., the (limiting) coexistence curve and diameter, 
even very close to T c . 



V. APPLICATIONS TO SIMULATION 

In this section we extend and illustrate the finite-size 
scaling analysis and the use of the special loci by esti- 
mating critical parameters for the hard-core square-well 
fluid [7] and the restricted primitive model electrolyte 
[19] on the basis of grand canonical Monte Carlo simula- 
tions. In particular, the Q-loci play an important role in 
determining the critical temperature and indicating the 
universality class of the models. Once the critical temper- 
ature is obtained, we may use the fc-loci and the Q-loci to 
estimate the critical density p c as already demonstrated 
in Sec. III.B: see Fig. 3. To estimate the universal cor- 
relation exponent v for the RPM, the critical isochore is 
then utilized. 



A. Estimation of T c for the hard-core square-well 
fluid 

The HCSW fluid is the simplest continuum model 
that exhibits realistic gas-liquid separation and critical- 
ity. Hard spheres of diameter a = a interact via an at- 
tractive square-well pair potential of depth e and range 
b = Xa. In the simulations discussed here [7], A is taken 
to be 1.5 which reasonably represents simple fluids such 
as argon, etc. Reduced temperature and density are de- 
fined, as usual, via T* = k^T/e and p* = pa 3 . 

As already observed, for systems with an axis of sym- 
metry, such as Ising ferromagnets and lattice gases, 
Binder [20,39] used the moment parameter Ul = (1 — 
1/3<3l) to estimate critical temperatures (and critical 
exponents) by evaluating the parameter as a function of 
T on the axis of symmetry, where, of course, the ordering 
field, h, vanishes identically for all L, and then locating 
self-intersections. However, asymmetric systems, such 
as continuum fluids where there is no obvious symme- 
try axis, pose a crucial question when one aims to apply 
the same idea: Where should one look? The best choice 
is, naturally, the locus of "symmetry" corresponding to 
the vanishing of the finite-size ordering field h(p, T, p; L). 
In practice, however, the mixing coefficients fci, ]2, and 
s 2 in the ordering field — see (1.4) and (1.8) — are not 
known for such systems so that it is difficult to determine 
the locus h = in, say, the (T, p) plane. 

Furthermore, suppose Ql is calculated along any fixed 
locus — such as the critical isochore or even, say, the lim- 
iting Q-locus Pq(T) — on which h does not vanish but, 
rather, remains nonzero for any L. The contributions to 
Q l from nonvanishing h may then be gauged by expand- 
ing the scaling function Q{xL,yL, • • •) in (4.40) about the 
critical point as 

Q{xl,Vl, ■■■) = Qc + QiXl + Qix\ + Q3VL 

+ Qavla + Q 5 y 2 L 5 + ■ ■ ■ , (5.1) 

where the linear terms yz,, 1/1,5, etc. vanish identically 
in view of the basic symmetry under yL^—VL, VL5^ 
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— UL5i etc. Evidently, any small uncertainties in the crit- 
ical parameters will be enhanced via the scaling combi- 
nation ul oc hL^l v when L increases. For example, if 
Sp c is an error in p c , the contribution to Ql will vary 
as y\ ~ 5p^.L 2 P/ v and hence diverge when L — > oo thus 
causing difficulties in extrapolating finite-size data. Ex- 
plicit calculations reveal the corresponding reduction in 
precision. 

Beyond this issue one finds, by explicit calculations 
for strongly asymmetric systems like the RPM, that the 
behavior of Ql{T) on the critical isochore, (p)l = p ci 
may not even be monotonic — as it is on the h = 
locus. This adds further uncertainty to interpreting the 
data. 

To overcome these obstacles, we consider the Q-loci, 
Pq{T] L), for a fixed L on which it was shown in Sec. IV. A 
that the scaling combination oc hL^l v actually decays 
as jiL~Pl v when L — > oo: sec (4.6). (Thus h vanishes like 
j2L~^ A+ ^/ v .) Hence, the Q-locus can be considered as 
an "optimal" choice for analyzing Q L and estimating T c . 
Notice, of course, that in a symmetric system the Q- 
locus reduces to h — (or, equivalently, to p = p c ). For 
certain other thermodynamic quantities one might find 
corresponding optimal loci, such as the /c-susceptibility- 
loci, etc. Here we examine Ql evaluated on the Q-loci 
for the HCSW fluid. (For application to the RPM: sec 
[19].) 

Generally, one must expect that Ql on a Q-locus starts 
near Q = | above T c (in the one-phase region); but, since 
the Q-loci in the two-phase region approach the diameter 
p(T) [see Fig. 6], Ql must then approach unity below T c 
[see (4.30)]. At T = T c the Q-loci approach the critical 
point so that Ql on a Q-locus must pass through the 
universal value Q c at some temperature, say T^{L), that 
approaches T c as L — > oo. These features are evident in 
the plots of Q on the Q-loci for the HCSW fluid shown 
in Fig. 10. Thus all the curves intersect one another 
near the Ising value Q c ~ 0.6236 (for periodic boundary 
conditions on a cube [21-23]) strongly confirming that 
the HCSW fluid belongs to the (d = 3)-dimensional Ising 
universality class. 

To obtain the asymptotic behavior of T®(L) for large 
L, wc solve the equation 

Ql(T;p q ) « Q(x L ,y L ,y L4 ,yL5,---)\Q = Qc (5.2) 

where the subscript Q notation denotes evaluation on the 
Q-locus. Substituting the expression (4.7) for jjl on the 
Q-loci and using (5.1), we may solve this equation to 
obtain 

x L = DltIL 1 ^ + • • • , 

« -U c L4 Q 4 /QiL e ^ - j 2 2 ^Q 3 /Qi£ 2/ ^, (5.3) 

where r was defined in (3.10) and Yq in (4.7), while the 
coefficients Qj in (5.1) could also be expressed in terms of 
the scaling-function expansion coefficients Yfi^. Finally, 
T®(L) is given by 
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FIG. 10. Plots of Ql(T; (p) L ) on the Q-loci, p Q (T; L), for 
the HCSW fluid providing estimates for T c and Q c - Classical, 
XY and Ising values of Q c are marked on the Q axis [21-23]. 
The system sizes match those in Fig. 6. 

#(L) = [T«(L)-T C ]/T C 

= -P X IL^I V - P 2 /L( 1+2 ^» + ■ ■ ■ , (5.4) 
Pi = QMJtQiDl, P2=jIY 2 Q 3 /Q 1 tD l . (5.5) 

Notice that for d = 3 Ising systems the leading exponents 
in (5.4) are (1 +6)/v ~ 2.41 and (1 + 2[3)/v ~ 2.62, 
the latter with an amplitude proportional to jf ; these 
large values explain the observed rapid convergence of 
the TQ(L). 

Figure 11 displays T®(L) versus for the HCSW 
fluid with the predicted Ising value ip = (1 + ff)/v ~ 2.41. 
The small value i? M = —,72/(1 — j 2 ) [see 1(3.41)] of about 
—0.04 discussed in Sec. III.B indicates that the ampli- 
tude Pi in (5.4) is negligible. Thus considering only the 
leading term is sensible. However, to allow for the vari- 
ous higher order corrections, the small shift parameter I* 
has been introduced. 

From this plot, we estimate the critical temperature 
for the hard-core square-well fluid to be 

T* ~ 1.2186 ±0.0003 (HCSW). (5.6) 

This value is about 0.06% higher than the estimate 
T* ~ 1.2179 ± 0.0003 of Orkoulas et at [7]. For the 
RPM, Luijten et al. [19] obtained a precision of ±0.04% 
in estimating T* by the same approach. 

It is worth stressing that in all these calculations (and 
those described above and below) it has been imperative 
to use extensive histogram reweighting procedures [44] 
in order to precisely determine intersections of loci, max- 
ima and minima, etc. It is clear that without sufficient 
precision and, indeed, accuracy in calculating finite-size 
properties, extrapolation procedures are doomed to fail- 
ure or, worse, seriously misleading estimates. 
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FIG. 11. Plots of T C Q (L) vs. (L* + l*)~* with 
ip = (1 + 0)/v = 2.41 to estimate T c * for the HCSW fluid. 



Unfortunately, the convergence of the estimates for Q c 
is not as rapid. We find 



Q% L (L) « Q c + {l + 9)Q i U c L jL e ^ 
+ (1 + 2(3)f 2 Q 3 Y*/LW», 



(5.9) 



where for Ising-type systems the exponents are 9/v ~ 
0.83 and 2/3/f ~ 1.04. This slower convergence may be 
the reason why the successive intersections seen in the 
inset in Fig. 10 suggest a limit some 1 or 2% higher than 
the established Ising value [21-23]. However, since no 
special efforts were originally made [7] to gather HCSW 
data optimal for evaluating Q and the Q-loci, one must 
also suspect the possibility of inadequate simulation ac- 
curacy. By contrast, the central unbiased estimate for Q c 
for the RPM (on which considerable effort was focussed) 
captured the Ising value precisely within uncertainties of 
only ±0.3% [19]. 



C. Estimating the correlation exponent 



B. Estimation of Q c 

There seems little serious doubt on the basis of Fig. 
10 (as well, of course, as on previous evidence [7]) that 
criticality in the HCSW fluid is of short-range Ising type. 
In other cases, however, one may well desire to estimate 
Q c , and hence resolve the universality class, in unbiased 
fashion. In that situation the successive intersections of 
plots of Ql on the Q-loci for increasing sequences of L 
values may be useful. Accordingly, let us define Tq L {L) 
and Qq L (L) as the intersections of a plot of Ql (T) on the 
Pq(L;T) locus with a plot of Ql-al{T) on the pq(L — 
AL;T) locus and ask for the asymptotic behavior as L 
increases at fixed, small AL. 

The analysis follows the lines of the previous section 
except that (5.2) is replaced by 

Q(xl,Vl, ■ ■ -)\q - Q(xl-al,Vl-al, ■ ■ -)\q »0. (5.7) 

For the temperature intersections we find 

\Tq L {L) — T c }/T c = OPi/ L^ 1+e ^ v 

+ 2/3P 2 /L (1+2 « /,y + • • • , (5.8) 

which, in leading order, is independent of AL. The coef- 
ficients Pi and Pi are the same as those defined in (5.5), 
that enter (5.4), namely the asymptotic result for t®(L), 
the intersections with Q c . However, the approach takes 
place from the opposite side, and since 8 ~ 0.52 and 
2/3 ~ 0.65, the amplitudes are smaller. For these reasons 
one might well prefer to use the successive intersections: 
however, a little reflection shows that they place greater 
demands on the precision and reliability of the simula- 
tions. 



Of basic importance and value in determining the uni- 
versality class of a model is the correlation length ex- 
ponent v. As already frequently stressed, this enters in 
finite-size systems via the combination L\t\ u which opens 
many routes to the estimation of v. For example, the 
scaling of Ql on the Q-locus should satisfy 

Ql(T; Pq{T; L)) - Q c AQ(tL 1/v ). (5.10) 

From this it follows that the derivatives, 
dQ L (T; p Q (T; L))/dT, evaluated at T c or at T®{L) or at 
Tq L (L), etc., will all, in leading order, diverge as L x l v . 
However, obtaining these derivatives accurately is a dif- 
ficult computational task. Furthermore, the corrections 
to the leading behavior are likely to be quite significant 
(owing, in particular, to the strongly nonlinear variation 
of AQ(x) which must saturate at constant values of order 
unity when x — > ±oo). 

To provide a robust method of estimating v from simu- 
lations above criticality — which are intrinsically easier to 
bring to equilibrium than simulations closer to or below 
T c — Orkoulas et al. [7] introduced various "estimator 
functions," yj(T,p). When evaluated in the thermody- 
namic limit on a critical locus, say £, that approached 
the critical point from above, these diverged as t — > 0; 
but in a finite system they exhibited rounded maxima 
above T c at temperatures Tj(L). For suitable loci, C, the 
Tj(L) must approach T c as L~ x l v . Then Orkoulas et al. 
considered unbiased exponent estimators, independent of 
the unknown (or known) value of T c . Specifically, for a 
pair yj and 34, they measured ATjk — Tj(L) — Tk(L) 
and computed sequences 

f ATjk {L + AL) 1 L 1 

AT jk (L) JAL^,' (5 - 11} 
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as L — > oo. By using estimates for the critical isochorc, 
Orkoulas et al. [7] estimated v for the HCSW fluid and 
confirmed its Ising-type character. They also checked 
that, within the available precision, the results for v were 
not sensitive to the estimate for p c . 

However, this method is relatively demanding in that 
the differences, Tj(L) — Tk(L), must be obtained to rel- 
atively high precision. For the RPM — which is much 
harder to simulate reliably than the HCSW fluid even 
above T c — this proved a stumbling block. In addition, 
while relative insensitivity to the estimate of p c could 
reasonably be expected, the very strong asymmetry and 
the likelihood of strong pressure mixing (since confirmed 
[25]) made the choice of critical locus more questionable. 
Would the critical isochore still be satisfactory? 

At issue in this latter question is that, as a result of 
pressure mixing, the estimator functions yj(T) pick up 
contributions varying with the fields h and p on the lo- 
cus £, say p — p c . This question is partly resolved by 
the analysis of I Sec. IV. D which shows that on the crit- 
ical isochore [and, by extension, on any locus behaving 
asymptotically as (p — p c ) « ct — > 0] one has h ~ |t| 1-0+7 
and p^ \t\ 2 ~ a . The associated correction exponents are 
sufficiently large (> 2) that they arc of little practical 
concern relative to the unavoidable leading correction- 
to-scaling terms varying as t e . In a finite system a dis- 
cussion along the lines leading to (5.9) (that invokes the 
analog of (4.7) for on the isochore) is appropriate; but, 
as in (5.9), the extra terms to be anticipated, varying as 
L~ 2/3 / u , are of higher order than the leading L~ B I V cor- 
rections. Nevertheless, it may be of value, as suggested 
in [9], to use as the locus ( a "theta locus" defined via 



MT) = Pc[0+(l-0)(Tc/T)], 



(5.12) 



where a most favorable value of d might be one chosen 
to approximate an optimal fc-locus or Q( fc )-locus. 

For the RPM a second problem arises which we explain 
here and then deal with explicitly. For completeness we 
recall that the restricted primitive model electrolyte con- 
sists of N = 2N + hard spheres of diameter a = a, of which 
N + carry a charge +qo and A_ (= N + ) a charge —qa. The 
pairwise Coulomb potential is iq^/Dr for two like/unlike 
charges at separation r. Appropriate reduced variables 
are 



T* = k B TDa/ql 



pa 



(5.13) 



Orkoulas et al. introduced twelve estimator functions 
y- (j = 1, • • • , 12) [7]. The simplest, y x = C v , was the 
constant volume heat capacity. But for the RPM this 
displays maxima fairly far below T c which, moreover, are 
not easy to locate precisely [45-47]. With 6 = 1/T* 
Orkoulas et al. defined y 2 — (<9CV/<90) p : this function 
has a local extremum, T 2 t "(L), above T c which varies fairly 
regularly as L increases: see Fig. 12. On the other hand, 
in the case of the RPM the functions 3^3, • • • ,3^6 prove 
to have maxima close to but below T c . The function 
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FIG. 12. The estimator function 3^2 (T) = (dCv /d&) with 
G = l/T* on the critical isochore (p* c ~ 0.079) of the RPM 
electrolyte (at a C = 5 discretization level [19]). The vertical 
line marks the estimated critical point at T c * ~ 0.05069 [19]. 



3^7, a modified susceptibility, displays no maxima on the 
critical isochore in the range 0.045 < T* < 0.070. The 
remaining functions 3^8 to 3^12 do display extrema above 
T c but their behavior is not very smooth for the accessible 
values of L*\ 

Accordingly, new estimator functions were sought. Af- 
ter some investigation two further acceptable functions 
were found, namely, 



2U/2 



d 2 (m 2 } 
dO 2 



6\l/6 



d 2 (m 6 ) 

de 2 



(5.14) 



where m= (N — (N))/V. The behavior of these func- 
tions resembles that shown for 3^2 (T) in Fig. 12 although 
for the same values of L the maxima lie further from T c : 
see K Figs. 4.15 and 4.16. 

Finally, we must accept that neither the quantity nor 
the quality of the obtainable RPM data suffice to im- 
plement the recipe (5.11). Instead, we accept the biased 
estimators 



A,- 



1 



Tj(L + AL) -T c 



ATj(L)-T c 



L*+l* 1 
AL* * z/ 



(5.15) 



which require a value for T c : that we take from the study 
of Q on the Q-loci as in Fig. 10 [19]. The shift parameter 
I* allows, as in Fig. 11, for higher order terms in the be- 
havior of the Tj(L). Extrapolation vs. 1/L, as illustrated 
in Fig. 13, yields 



v = 0.63 ±0.03 (RPM). 



(5.16) 



This value (previously reported but not justified [19]) 
supports the conclusion that, despite the infinite range 
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FIG. 13. Plots of the estimators (5.14) for the exponent 
1/u for the RPM using y 2 (open circles) and y' 6 (crosses) 
with I* — 2, 0, and —2 from the top downwards, and y' 4 (solid 
circles) with I* = 5, 3, 1, —1, and —3: see text and (5.14). 

of the ionic forces underlying the model, it behaves, as re- 
gards phase separation and criticality, like a short-range 
Ising-type system. 

VI. FULL SCALING IN THE CANONICAL 
ENSEMBLE 

In the thermodynamic limit for regular systems 
there is a full equivalence between the different en- 
sembles. Consequently a "canonical description" in 
terms of the Hclmholtz free energy density f{p, T) — 
liniL^oo F N (V,T)/V with p = lim L ^ 00 (7V/F), is as valid 
and provides the same information as the grand canon- 
ical viewpoint based on p(T, p) that we have so far 
adopted. Similarly, as observed in the Introduction, in 
leading order the canonical scaling form (1.13), which 
invokes the scaled combination z oc m/\t\^, is equiv- 
alent to the grand canonical form (2.4) which entails 
y oc h/\t\ A oc [p — p a (T)]/\t\ A . However, in higher orders 
the necessity for field mixing via (1.1)-(1.4) complicates 
matters. Specifically, whereas the full scaling fields t, p 
and p are generally nonsingular functions of the under- 
lying scaling fields, t, p, and p (unless renormalization 
group "resonances" arise [27]), this is no longer the case 
for the canonical variables fa, t and /. Here we derive 
some of these complications that arise canonically, first 
in the thermodynamic limit in the presence of pressure 
mixing, then in finite systems. In the latter case we wish, 
in particular, to understand the asymptotics of the finite- 
size, classical-type critical points that may be identified 
in canonical simulations: see, e.g., [7,19]. 



A. Thermodynamic limit 

By standard thermodynamics for infinite systems the 
Helmholtz free energy density is given by 

f(p,T)=pp-p, (6.1) 

where p and p are understood to be re-expressed in terms 
of the density via p = (dp/dp)T- It is straightforward to 
introduce the reduced variables p, p and 

(6.2) 

via (1.1) and (2.12), and convenient to recall (2.17), for 
ei and and, further, to write 

e = e 1 " 1 , e 2 =j 1 +j 2 l 1 , e 4 = fci+j 2 fco- (6.3) 

Now we must address the choice of general canonical 
scaling variables. We wish, first, to allow for the leading 
correction-to-scaling terms which are expressed in terms 
of t both for infinite and finite systems in (2.2)-(2.4). 
Accordingly, it seems appropriate to adopt t also canon- 
ically, although it will need to be re-expressed in terms 
of ra in place of p. 

Similarly, it seems clear that the general scaling field fn 
should be chosen conjugate to the general ordering field 
h. Thus we adopt 

m = (dp/dh) h (6.4) 

which is identical to the scaling density p that was intro- 
duced in (2.10) along with the scaling entropy s. 

With these variables in hand we can rewrite (6.1) as 

f(p, T) ee f(p, T)/ Pc k B T c = Up, T) + f(p, T), (6.5) 

where the nonsingular background term may be ex- 
panded as 

fo{p, T) = f c + p c p - fc t + e e 4/ 5t H , (6.6) 

with f c = {pcPc-Pc)/p c kBT c and p c = p c /k B T c . On the 
other hand, the singular contribution becomes 

f(p, T) = rhh — p + j 2 mp - e e 3 sh - j 2 eoe 3 sp H , 

(6.7) 

in which the presence of the coefficient j 2 makes clear 
how pressure mixing enters. 

Our aim now is to express / in terms of the general 
canonical scaling combinations 

z = m/B\tf, yi = U 4 \t\ e , y 5 = U 5 \tfs, (6.8) 

where B = QU: see (2.4), (2.5) and accompanying text. 
To that end, from (2.4) and (6.4) we first obtain 

z = W^(y;y 4 ,y 5 ,---) with y = Uh/\t\ A , (6.9) 



p = e (h+j 2 p + e 4 t H ), 



19 



for t ^ 0, where W'±{y; ■ ■ •) = dW±/dy. Inverting this 
expression yields 

h=U- 1 \t\ A F£(z;y 4 ,y 5l ---), (6.10) 

where the scaling functions F±(z) are the inverses of the 
W±(y). From (2.4), we hence find 

p = Q\t\ 2 - a F£(z;y 4 ,y 5r --), (6.11) 

in which the new scaling functions are defined by 

F£(z;y4,v B ,---) = W ± (F£{z;y i ,---);vi,---). (6.12) 

Then, rearranging (1.2)-(1.4) and substituting yields 
the canonical thermal scaling field as 

t = rt- e e 3 h - e Q e 2 p H , 

= Tt-(e e 3 /U)\t\ A F£(z;y 4 ,---) 

- e e 3 Q\t\ 2 - a Fl(z; y 4 , .••) + ••• , (6.13) 

where r = 1 — e (fco e 2 + k dso defined in (3.10). 

When the mixing coefficients, l\, j\, and j 2 all vanish t 
reduces to rt. Notice, however, in contrast to the grand 
canonical formulation, that for nonzero l\ or j\ the scal- 
ing fields t is now a singular function of t with leading 
nonlinear contributions varying as \t\ 2 ~ a ~P and \t\ 2 ~ a 
( in place of t 2 , etc.). By the same token, corrections 
proportional to fh 2 ~ m 2 , arising from the expansion of 
F±(z) and F±(z), will carry the singular factors |i| 7- ^ 
and i| 7 ; moreover, the former actually dominates the 
nominally leading term linear in t. 

For the general canonical order variable, fh, we find 
from 1(2.18) 

fh = e a m + e e 3 s + (j 2 + jifci)(e§/r)m 2 H , 

= e Q m + e Q e 3 Q\t\ 1 ~ a Fl(z; y 4 , ■ ■ •) 

+ (j2+nki)(e 2 /T)m 2 + -.., (6.14) 

where, from (2.10) for s, we find 

Fi(z; • • •) = (2 - a)*£(z; •••)-(/? + 7)^(^5 ' ■ ■)■ 

(6.15) 

Evidently, fh also entails singular terms which, indeed, 
introduce 1 1 as a leading correction unless e 3 = h+ji 
vanishes. 

Finally, f(p, T), the singular part of the Hclmholtz free 
energy, can be expressed as a sum of a scaling piece, 
which simply extends the original leading form (1.13), 
plus a series of nonscaling, singular but higher order cor- 
rections arising from field mixing. If we define the scaling 
functions 

X ± (z;y4,---)=FZ-zF£, X p ± = zF£, (6.16) 
X^{z- • ■ •) = ± FUz; ■ ■ -)F£(z; • • •), X' ± = ± F^F%, 

the explicit result, recalling (6.8), is 



f(p,T) = -Q\t\ 2 - a [X ± (z;y 4 ,y 5 ,---) 
- j2 QU\tfX p ± (z;---) 

+ (eoe 3 /U)\t\ 1 - a -Px£(z;---) 

+ j2e a e 3 Q\t\ 1 - a X s ± (z; ■■■) + ■■■]. (6.17) 

Evidently the most singular nonscaling correction is of 
relative order \t\@ and arises only from the pressure mix- 
ing coefficient j 2 that induces a Yang- Yang anomaly. In 
as far as this and the other nonscaling corrections are 
of higher order (in powers of t) than the scaling term, 
they might be regarded as part of a "singular background 
piece", say, fo s (p,T). But the singular nature of the 
canonical scaling fields t and fh cannot be so readily 
sidestepped! 

In contemplating these results one may speculate that 
there might exist better choices of the canonical scaling 
fields, fh and t, that would ameliorate the singular mix- 
ing terms in (6.13) and (6.14) and/or absorb some or 
all of the nonscaling corrections in (6.17); however, this 
seems unlikely to us. Indeed, it is worth recalling that 
even the concept of a "nonsingular background" encoun- 
ters dangers near criticality in a canonical or Hclmholtz 
formulation. Thus in a symmetric system near T c with 
m = p — p c one might reasonably expect the background 
to have the power series expansion 

f (p, T) = f c + + f 2 m 2 + f h2 tm 2 + f 4 m 4 + ■■■. 

(6.18) 

But since the inverse susceptibility x _1 (^) is given by 
(d 2 f/dm 2 )T, the susceptibility itself cannot diverge at 
T c unless f 2 vanishes identically. Similarly, if fi t2 and f 4 
do not also vanish one would have 7 < 1 and 8 < 3, both 
of which inequalities contradict exact theory and precise 
experimentation! These observations point, of course, 
to the fundamental character of a grand canonical or, 
better, a full field formulation in terms of p, p,, and T. 

B. Finite-size canonical criticality 

To extend our canonical scaling description to finite 
systems we may follow Sec. II. First, in the set of scaled 
variables (2.3), we replace j/x, = UhhL^I 1 ' by 

z L = B L fhL () l v . (6.19) 

Then, in addition to a nonsingular background free en- 
ergy fo(p,T;L), we may anticipate a singular part, cor- 
responding to (6.17), of the form 

f s (p, T; L) = L- (2 -^^[X (x L ,z L ; y L4 , y L5 , ■ ■ •) 
+ .hL^ l3/ly X 1 (x L} z L ;y L4 , ■ ■ •) 
+e e 3 L^-^^X 2 {x L ,z L ----) 
+ h^e 3 L^l v X 3 {x L , z L - •••) + •••] • (6.20) 
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Note that the new finite-size scaling functions Xq and 
X 3 should be symmetric under zl^-zl, J/ls — 2/5L, 
etc., while X\ and X2 are antisymmetric. In the absence 
of field mixing we recover the obvious finite-size general- 
ization of the scaling form (1.13). However, the pressure 
mixing coefficient j 2 generates a nonscaling correction 
that vanishes as L~P/ V and is antisymmetric in zl, yL5, 
etc. The coefficient l\, that mixes the chemical poten- 
tial into the thermal field t, produces an antisymmetric 
correction vanishing as L~ { - 1 ~ a ~P s >l v . 

As in (2.6)-(2.8) we expect that the scaling functions, 
Xj(xl, Zl', Vl4, Vl5, • • can be expanded generally in 
powers of the irrelevant variables, ula, UL5, etc., and, 
also for finite L near criticality, in powers of xl and zl, 
with coefficients Xf kl as in (2.8). [See also K(4.182)- 
(4.184).] There is, in fact, a concealed subtlety here: 
specifically, the particle number N is an integer so that 
the density p (and m) are intrinsically discrete variables 
in a finite system. Away from criticality the free energy 
surely approaches an analytic function of p when L — > 00; 
but the degree to which a corresponding smoothness may 
be assumed in a finite system close to criticality is not 
obvious. [Incidentally, the corresponding issue can be 
raised in connection with the two-Gaussian description 
of the distribution Pl(p', p, T) in (4.13).] However, in the 
absence of concrete evidence to the contrary, the assump- 
tion that the finite-size canonical free energy, f(p,T;L), 
may be treated as an analytic function through (p c ,T c ) 
seems highly plausible if used, as here, to determine lead- 
ing asymptotic behavior when L — > 00. 

Now simulations of simple fluid systems reveal that 
as, a function of density, f(p, T; L) exhibits two peaks 
for T < T c that correspond to the separation of the 
two phases. One may then define a finite-size canoni- 
cal critical point, (p®(L), T C °(L)), as a point where these 
two peaks merge. By virtue of the analytic behavior of 
f(p,T;L), such canonical critical points must, in gen- 
eral, be classical in character. However, they will — at 
least in simple cases — approach the bulk critical point 
(p c , T c ), whether or not the critical behavior remains clas- 
sical in the thermodynamic limit. In principle, extrapo- 
lating such canonical critical points may help locate the 
limiting critical point; in practice, however, this has so 
far proved of limited usefulness [7,19]: see the numeri- 
cal behavior revealed in Fig. 3 of [7] and Fig. 1 of [19]. 
Nevertheless, it is of interest to elucidate the asymptotic 
behavior, especially of p®{L). 

The conditions determining a classical critical point 
reduce to 

(df s /dm) T = 0, (d 2 f s /dm 2 ) T = 0. (6.21) 
On expanding the scaling functions in (6.20) these yield 

= 2X% 2 z L + ]2 X% 1 L-PI' J + 2X { Q %U c Li L- e /»z L 

+ e e 3 Xl 01 L^^ + ---, (6.22) 
= 2Y ° 02 + 2Xg <12 x L + 6j 2 Xl 03 L-^»z L 



+ 2xl ) %Ul i L- e ' v + ■■■. (6.23) 

Solving these equations for xl and zl and using (6.13) 
and (6.14) for t and m finally yields the critical temper- 
ature as 

t c (L) = lT?(L)-T c ]/T c 

= c x L- x l v [\ + c 2 L- e ' v + 32C 3 L- 2S) I V + ■■■} (6.24) 

and the canonical critical density as 

p° c (L) = Pc [l + jihL-W" + feL^ 1 -")/" 

+ b 3 L-W +e V v + ■■■}, (6.25) 

where the leading amplitudes are given by 

ci = -K02/ DtX° A2 and h = -e 1 X 1 ° 01 /2X ° 02 . 

(6.26) 

It is instructive to learn that the asymptotic behavior of 
Pc(L) has the same form as exhibited by the fc-loci and 
the Q-loci evaluated at T = T c : see (3.13) and (4.8). 
The data for the RPM, however, suggest that the two 
leading corrections in (6.25) compete rather strongly so 
thatp°(L) appears to approach p c nonmonotonically [19]. 

VII. CONCLUSION 

In this article we have extended to finite systems the 
"complete" scaling theory developed in Part I [10] for 
critical behavior in the thermodynamic limit that incor- 
porates pressure mixing in the scaling fields as well as 
the irrelevant corrections to scaling. The basic theory is 
set out in Sec. II in a grand canonical or (p, p, T) formu- 
lation: sec (2.1), (2.3), and (1.1)-(1.4). The possibility of 
finite-size corrections in the scaling fields p, h, and i [see 
(1.8)] has been reviewed briefly in Sec. II. B and, in Sec. 

II. D, a fairly direct route to detecting such a dependence 
- by studying numerically p(T c , p c ;L) — is proposed. 

Section III applied the theory to elucidate the near- 
critical behavior of the fc-loci, defined in the (p, T) plane 
by the isothermal maxima of the modified susceptibili- 
ties x{T,p)/p k - see Fig. 1. The usefulness of the fc-loci 
in estimating the critical density, p c , via simulations is 
demonstrated for the hard-core square-well fluid in Sec. 

III. B and Figs. 2 and 3. It also transpires that the value 
of k which yields a locus that approaches the critical 
point "mostly directly" provides a reasonable estimate 
of the Yang- Yang ratio i? M [8-10] that, in turn, provides 
the most direct measure of the degree to which pressure 
enters the ordering field h. In this way Fig. 1(b) provides 
rather clear evidence of a significant ratio, ~ 0.26, in 
the restricted primitive model electrolyte: see Sec. III.B. 

The behavior of the basic moment ratio Ql(T; (p)), as 
defined (following Binder [20]) in (1.9), is the topic of Sec. 
IV: see Figs. 4 and 5. In particular, the associated Q-loci 
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(and QW-loci) are determined in Sec. IV. A (and IV. B): 
sec (4.8) [and (4.10)] and Fig. 6. Of especial interest 
is the behavior of Ql(T; (p)) below T c , within, up to, 
and beyond the boundaries, p+(T) and p~(T), of the 
two-phase region: see Figs. 7-9. For fixed T < T c and 
large enough system sizes, L, exact nontrivial results have 
been found, as shown in Sec. IV. C and D: in particular, 
the study of the minima in Qr,{T,(p)) [see Fig. 9 and 
(4.38)] lays the foundation for a precise method [25] of 
estimating [p + {T) — p~(T)] and the coexistence curve 
diameter p{T) at higher temperatures very close to T c . 

Of remarkable value for estimating T c for asymmetric 
fluid models is the behavior of Ql evaluated on the corre- 
sponding Q-loci: see Fig. 10 and the asymptotic expres- 
sion (5.4) and corresponding plots in Fig. 11. Likewise, 
the estimation of the critical value Q c = Qoo(T c ; p c ), de- 
scribed in Sec. V.B, is important for determining the uni- 
versality class of criticality. Finally, in Sec. V.C and Figs. 
12 and 13, the estimation of the critical exponent v for 
the highly asymmetric restricted primitive model elec- 
trolyte has been described (confirming Ising character). 

The issue of a canonical or (p, T) formulation of crit- 
icality with corrections to scaling and pressure mixing 
is taken up in Sec. VI. The basic expression, (6.17), for 
the singular part of the Hclmholtz free energy is intrin- 
sically more complex than the (p, p, T) scaling formula- 
tion, entailing an infinite series of "improperly scaling" 
corrections. This formulation provides a basis for deter- 
mining the asymptotics of the canonical critical points 
(of classical character) that can be observed in (N, V, T) 
simulations: see (6.24) and (6.25). 

In summary, we believe that the theory developed here 
and the applications illustrated constitute a solid foun- 
dation for future computational studies of criticality that 
employ systems of finite size. 
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